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We apply the DMRG method to the Pariser-Parr-Pople hamiltonian and investigate the onset 
of dimerization. We deduce the parameters of the hopping term and the contribution of the a 
bonds from ab initio calculations on ethylene. Denoting by Rij the C — C distances, we perform a 
variational optimization of the dimerization S = — Ri-\^)/2 and of the average bond length 

Ro for chains up to N — 50 sites. The critical value of N at which the transition occurs is found to 
be between N = 14 and N = 18 for the present model. The asymptotic values for large N for Rq 
and S are given by 1.408(3)A and 0.036(0)A. 



I. INTRODUCTION 

The quantum chemical treatment of large polymeric 
systems is extremely demanding from the computational 
point of view. When the electron correlation is taken into 
account, traditional correlated quantum chemical meth- 
ods grow too rapidly with the number of electrons (n 5 
or worse) to be considered as practical tools to study 
large systems. Density Functional Theory (DFT) is a 
powerful tool to cope with this kind of problems, and a 
number of efficient codes are available nowadays. How- 
ever, the local density approximation (even with gradient 
corrections) can only treat the local part of the electron- 
electron correlations (within an atom or bond) and does 
not treat the interatomic correlations with sufficient ac- 
curacy. Moreover, no information on the wave function of 
the system is provided by this method. In this paper we 
focus our attention on a different method, i.e. the Density 
Matrix Renormalization Group (DMRG) formalism, pro- 
posed by White! This approach in principle can also be 
developed to an efficient tool to compute correlated ener- 
gies and wave] 'unctions of the ground or excited states of 
large systems. DMRG is an approximation method that 
uses a density matrix to identify the most relevant states 
of a full configuration interaction (FCI) expansion. The 
localized orbital set is divided into two adjacent parts 
called "system" and "environment" . The DMRG method 
allows one to determine in a precise mathematical way 
the most probable states of the " system" in the presence 
of the "environment". These states are retained, and the 
others are dropped. The calculation starts with two small 
fragments, whose size is increased in the course of the 
calculation, until the whole system reaches the wanted 
dimension. We refer the reader to Refs. [ and [ ||] 
for a detailed description of the method. 

DMRG has been applied to a wide class of systems, in- 
cluding strongly-correlated electrons, antifcrromagnetic 



Heisenberg chainso,jphononstH3, defined on one and two 
dimensional latticeso. Originally DMRG was designed 
to treat nearest-neighbor interactions in one dimension, 
and in these conditions the method is very efficient in 
reducing the computational effort. 

The Pariser-Parr-Pople (PPP) hamiltonian of fippju- 
gated polyenesS was recently studied with DMRGcffl. In 
Ref. [ |^| the dimerized problem.. was considered with a 
symmetrized version of DMRGEl that allows the selec- 
tion of excited states of given symmetry. In Ref. [ ||] the 
ground state at constant bond length was examined, and 
the results were in full agreement with FCI calculations, 
in all cases where these latter were availableBB. The 
PPP hamiltonian has a long history since in spite of its 
simplicity it provides at least a qualitative explanation 
of the essential physics of polyacetileneEao. We decided 
to consider the dimerization of the Pariser-Parr-Pople 
model of conjugated polyenes. Systematic studies of the 
bond alternation in conjugated polyenes were carried out 
by Paldus and coworkersE£lt3. These authors used vari- 
ous kinds of approximations (self-consistent field, alter- 
nant molecular orbitals, coupled cluster, valence bond) 
since the FCI treatment ceases to be applicable before 
the onset of the bond alternation. It seems therefore 
interesting to compare those results with DMRG calcu- 
lations on the same hamiltonian. Recently a Density 
Functional (DF,I|) computation on large annulenes has 
been performedtll, so a comparison of our DMRG results 
obtained with a model hamiltonian with a realistic all- 
electron calculation can give useful information on the 
properties of the PPP model. DMRGr-has already been 
applied to the study of dimerizationllaEll of polyacetilene 
chains using various model hamiltonians. In Ref. [ [l8| , 
the authors considered a generalization of the Hubbard 
model to describe the valence 7r electrons and approxi- 
mated the a contributions by a sum of pair interactions. 
Their generalization of the Hubbard model consisted of 
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allowing for a distance-dependent hopping term, while 
the Coulomb on-site repulsion was taken to be indepen- 
dent of the distance between C — C atoms. They de- 
termined the parameters of the model from accurate ab 
initio CI calculations on the ground state and the first 
excited state of the ethylene, molecule, available as a func- 
tion of the C — C distanceca. The dimerization 8 and the 
average bond-length R were determined minimizing the 
energy. 

The aim of this paper is to exploit the same method 
for the determination of the parameters as a function of 
the C — C distance, and to improve the treatment of the 
correlation by substituting the crude Hubbard approx- 
imation with the more accurate PPP treatment of the 
Coulomb repulsion. For this purpose a generalization of 
the PPP hamiltonian is used where both the hopping in- 
tegral and the potential depend on the geometry of the 
system. In this way, the values of the average C — C dis- 
tance and of the bond-alternation can be determined by 
a variational calculation as a function of TV. Furthermore 
different geometries can be easily considered, which was 
not the case for the Hubbard model. 



II. DEFINITION OF THE MODEL 

To describe dimerization in trans-polyacetilene chains 
we introduce the following distance dependent hamilto- 
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H = H„ 



Hi, 



(2.1) 



The a contribution is approximated by a sum over all 
bonds and depends only on the bond lengths: 



H a — E E a (Rij) 



(2.2) 



<hj> 



(Rij is the distance between two carbon atoms located 
at sites i and j and < ij > denotes summation restricted 
to nearest neighbors). The ir contribution is represented 
by the Pariser-Parr-Poplc hamiltonian, with distance de- 
pendent hopping: 




FIG. 1. Structure I: all C atoms are on a ring. 
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* = E & A + \ E T« (fk - 1) (»j - 1) (2-3) 



<i]> 



i,j=0 



Here are the generators of the unitary group 
summed over spin, and hi = En is the occupation num- 
ber of the site %\ /3y, 7y are parameters of the model. 
For the Coukimb repulsion we use the Mataga-Nishimoto 
prescription^: 
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7« 



7o 



Ri. 



(a. u.) 



(2.4) 



where 70 = 10.84 eV. 

Following Ref. [ |is}j we determine (3ij and H a by impos- 
ing that the singlet triplet splitting calculated ab initio 
for the ethylene molecule is correctly reproduced. In Ref. 
[ |22| the singlet-triplet splitting and the sigma energy are 
given as polynomial expansions in the bond length: 

5 

g(R) = ( 3 E(R) - 1 J B( J R))/2= £ ai R l (2.5) 

{=1 

5 

E a (R) =Y, a 'i Rl ( 2 - 6 ) 

(the coefficients a/ and a[ are given in Table III of Ref. [ 
p2]). To make contact with our hamiltonian, we calculate 
the singlet triplet splitting applying H„ to the ethylene 
molecule. Let us denote the singlet and triplet lowest 
energy states as: 
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(2.7) 
(2.8) 



It is easy to show that 3 E = (t\H v \t) = while the singlet 
energy is: 



i E(R) = {s\H„\s) = ±[( 7o - 7h2 (R)) 
-y/i-yo -7i,2(i?)) 2 + 16/? 1 2 2 (i?) 



(2.9) 



where the 7y are given in Eq. 2.4. This is to be equated 
to —2g(R) and solved to give the dependence of (3 = 
on the bond length: 
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(2.10) 



Using the Mataga Nishimoto prescription for 70, we ob- 
tain for example (3 = -2.936160eU when R = 1.40l. 
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FIG. 2. Structure II: all C atoms are on the surface of a 
cylinder, all bond angles are equal to 120°. 



III. DMRG ALGORITHM AND BOUNDARY 
CONDITIONS 

We are now ready to implement the DMRG procedure. 
At every stage the actual sites are decomposed into sub- 
blocks forming a " superblock" A • • B (where • denotes 
a single site). Initially A and B are single site blocks. In 
order to enlarge the blocks, the reduced density matrix of 
the "system" A • is computed and its most representative 
states are retained to define a basis of the new block A 
for the next iteration. The blocks B are obtained from 
the blocks A through reflection. 

Once the total size of the whole system is reached, 
the algorithm may be continued with an asymmetric de- 
composition A • • B with A and B different from one 
another. We can increase A and decrease B and vicev- 
ersa, optimizing the states that describe the blocks. This 
part is called finite-size algorithm, opposed to the first 
part called infinite-size algorithm (see Refs. [ for 
details). 

Generally in DMRG calculations open boundary con- 
ditions are imposed because the decay of the density ma- 
trix eigenvalues is faster with respect to the case of peri- 
odic boundary conditions and therefore the convergence 
is more accurate. We have tested that this difference is 
still present even with a long ranged potential, as in our 
case. When the system is finite the properties of the sys- 
tem depend on the choice of the boundary conditions. We 
put our system on a ring (structure I shown in Fig. [I]), 
since it is a common experience in finite-size calculations 
that the absence of boundaries ensures a nice (smooth) 
scaling of the energy with the system size. Let us denote 
by r c the circle radius. We want to consider only closed 



shell systems, since the behaviour of the .closed shell en- 
ergies with respect to N is quite smoothEa (the same is 
also true for low lying excited states that have an open 
shell structure and the same symmetry). 

To obtain a closed shell state, we use periodic bound- 
ary conditions (PBC) for N = Av + 2 and antiperiodic 
boundary conditions (ABC) for N = Av (denoting by 
Ci the annihilation operator corresponding to the site i, 
ABC amount to put cn — —cq). In fact the allowed 
pseudomomenta are: 



%2n 



PBC 

|(2n + l) ABC 



We denote by aoi and «i2 the two angles seen from 
the center corresponding to the double and single bond 
respectively, by the angle between sites i and j, by 
Rq the average bond-distance and by 6 the dimerization: 



R 



Ri 



S = 

From the equations: 

Rq — S = 2r r _ sin 
Ro + S 



aoi + a>i2 



2r c sin 

2tt 
~N 



m l — 

. /ai2 

m 

V 2 



(3.2) 

(3.3) 
(3.4) 
(3.5) 



it is easy to deduce that the circumference radius is given 
by: 



i? 2 cos 2 ( . )+J2sin 2 ( . ) 



sin(f) 



(3.6) 



FIG. 3. Structure III: annulene geometry, proper of 
6(2n + 1) systems, here shown for N = 18. All bond angles 
are equal to 120°. 
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Other geometries are also possible. In particular, one 
can fix the bond angles to 120°. This value is close to the 
experimental bond angle of the sp 2 hybridized carbon in 
linear polyenes. Two different geometrical arrangements 
of the atoms with bond angles of 120° are shown in Fig. ^ 
and|3[ Imposing periodicity on a linear chain implies that 
atoms are disposed on the surface of a cylinder (structure 
II shown in Fig. |^). The cylindrical geometry can easily 
be implemented within the above DMRG scheme. The 
annulene D^h geometry, proper of 6(2n + 1) systemsEZl 
(structure III shown in Fig. |3|), requires a slight mod- 
ification of the DMRG superblock arrangement, which 
must be changed to uABu due to the lower symmetry 
(translation invariance is lost). We note that the sigma 
part of the PPP hamiltonian does not depend on bond 
angles. Therefore the energy difference between the var- 
ious geometries is purely due to the distance difference 
between second and further neighbours in the molecule. 

IV. ACCURACY TEST 

There are two ways of improving the accuracy of a 
DMRG calculation: increasing the number m of states 
kept to describe A and performing one or more iterations 
of the finite-size algorithm. For comparison, we compute 
two sets of values, obtained respectively with the infinite- 
size algorithm and m — 256 states (DMRG(l)) and with 
the finite-size algorithm and m = 512 states in the last it- 
eration (DMRG (2)). Even when we limit ourselves to the 
infinite-size algorithm we must use from the very begin- 
ning the hamiltonian coefficients that describe the final 
system, because the distances depend on the final size. 

We test both methods against FCI calculations for 
TV = 14. As Table | shows, the absolute error of 
DMRG(l) is of order 1 x 1(T 2 eV, while the absolute 
error for DMRG(2) is of order 2 x 1(T 4 eV. 

TABLE I. Energy results for N = 14: the tv energies calcu- 
lated with DMRG are compared to FCI values. In DMRG(l) 
we use the infinite-size algorithm and keep 256 states. In 
DMRG(2) we use the finite-size algorithm and keep 512 states 
in the third (last) iteration. Ro and 8 are in A, energies in 



eV. 



Ro 


8 


E n [DMRG(l)] 


E n [DMRG(2)} 


E4FCI] 


1.390 


0.0000 


-35.678793 


-35.685545 


-35.685710 




0.0100 


-35.718946 


-35.725648 


-35.725808 




0.0200 


-35.835913 


-35.842554 


-35.842695 


1.400 


0.0000 


-34.873826 


-34.880260 


-34.880426 




0.0100 


-34.913695 


-34.920161 


-34.920230 




0.0200 


-35.029901 


-35.036097 


-35.036154 


1.410 


0.0000 


-34.083595 


-34.089740 


-34.089900 




0.0100 


-34.123179 


-34.129256 


-34.129409 




0.0200 


-34.238365 


-34.244233 


-34.244367 



Nonetheless, dimerization and average bond length are 
in substantial agreement: a polynomial fit of the data 
(even in S) gives R = 1.404141 for DMRG(l) , R = 
1.404101 for DMRG(2) and FCI, and S vanishing in all 
three cases. Anyway we do not learn much from the 
N = 14 case since the difficult task is the determination 
of 6. We note that the total energy is not very sensitive 
with respect to variations of 5. 

For N = 16, S is non-vanishing. In this case we test 
the accuracy of the method by comparing energies ob- 
tained exchanging double and single bonds. The hamil- 
tonian with (anti)periodic boundary conditions is invari- 
ant with respect to the exchange of double and single 
bonds. Our arrangement used in DMRG is not invariant: 
the bonds connecting the blocks are always asymmetric 
with respect to exchange of single and double bonds and 
also the bonds internal to blocks A and B are asymmet- 
ric for certain numbers of sites. So the symmetry will 
only be recovered if convergence is accurately achieved. 
Let En(+5) (E n (—S)) denote the 7r-energy obtained with 
i? 0j i a double (single) bond. Table || shows these sets of 
energies for the case N = 16. 

We find that DMRG(l) does not provide a very ac- 
curate determination of 5: minimizing over the set of 
energies E(+S) we obtain the values Ro = 1.405931, 
5 = 0.02361, while minimizing over the set of ener- 
gies E(—S) we get a quite different minimum Rq = 
1.405991, S = 0.02401. DMRG(2) improves the re- 
sult considerably: from the set of energies E(+S) we get 
R a = 1.405771, 5 = 0.0227(6)1, from E(-5) we get 
R = 1.405771, 5 = 0.0227(2)A Thus we see that the 
finite-size algorithm proves to be very useful. In order 
to spare computing time in this case, it is worth imple- 
menting a transformation of the wavefunction that gives 
a guess-ifor starting the next iteration, as suggested by 
WhiteEi 

TABLE II. Energy results for N = 16: comparison of 
DMRG(l) and DMRG(2) 7r-energies E v (+8) (E n (-S)) ob- 
tained by choosing i?o,i a double (single) bond respectively. 
E n (+8) and E^(— 8) are very close for the case DMRG(2). Ro 
and 8 are in A, energies in eV. 

DMRG(l) DMRG(2) 

Ro 8 Evj-8) E^j+S) Ek(-8) E n (+S) 

1.400 0.020 -39.908840 -39.908750 -39.920105 -39.920102 
0.025 -40.014971 -40.014343 -40.024754 -40.024635 
0.030 -40.138842 -40.138400 -40.147548 -40.147519 

1.405 0.020 -39.456517 -39.456962 -39.467522 -39.467413 
0.025 -39.562115 -39.561796 -39.571654 -39.571722 
0.030 -39.685307 -39.685028 -39.693748 -39.693736 

1.410 0.020 -39.008484 -39.008540 -39.019108 -39.019058 
0.025 -39.113436 -39.112946 -39.122670 -39.122689 
0.030 -39.235946 -39.235557 -39.244184 -39.244095 
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FIG. 4. Total energies at fixed Ro = 1AA. Results are 
obtained with DMRG(2). S vanishes for N < 14. 



V. RESULTS AND DISCUSSION 

In Fig. ^ we show the behaviour of the energy versus 
S for fixed R = 1.40l, using DMRG(2). In contrast 
with the case of open boundary conditions used in Ref. 
[ [l8), we find \s=o ~ 0, as expected. We observe the 
transition from a non-dimerized minimum for N < 14 to 
a dimerized minimum for N > 18. 

We perform two-dimensional minimization (energy 
versus 5 and R ) with both DMRG(l) and DMRG(2) 
as shown in Fig. and|^. We find a discontinuity in the 
slope of the curve of Rq versus N in correspondence of the 
transition from a non-dimerized minimum to a dimerized 
one (N > 14). 
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FIG. 5. Average bond length and dimerization, as de- 
termined from minimization of total energies with DMRG(l) 
(m=256 states) and DMRG(2) (m=512 states). Data denoted 
with a star indicate values of the average bond length obtained 
from minimization with the constraint 5 = 0. 



In fact it can be seen from Fig. || that the curve of 
i?o values obtained fixing § = shows no discontinu- 
ity. The asymptotic values for R$ and S are approxi- 
mately given by 1.408(3)1 and 0.036(0)1. The critical 
value N = 14 was also found by Paldus and coworkers 
in his VB calculationEj as well as in the previous studies 
of the series on dimerization of PPP cyclic polyenescJ. 
The hamiltonian used has the same two-body part, but 
we use a different prescription for the dependence of the 
hopping integral /3 and the sigma energy E t7 (R i j) .upon 



the annulenes CnHn, with N < 66 give a critical value 
N =j—30, confirming and extending previous results of 
MP2E3 and DFT calculations^. A direct comparison of 
the results is rather problematic due to the difference 
between the models. In our case we have a simplified 
hamiltonian treated at a high level of approximation by 
the DMRG method. On the other hand, the DFT cal- 
culation relics upon a more realistic all electron ab-initio 
hamiltonian, with a poorer treatment of the long ranged 
electron correlations. 
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FIG. 6. Single and double bond lengths as obtained from 
energy minimization with both DMRG(l) and DMRG(2). 

The location of the minimum of the energy as a func- 
tion of R and S is the result of a balance of the E„ 
and E v energies, as can be seen from Fig. [z| and one 
might wonder whether a small variation of one of the 
two opposite terms could significantly change the results. 
The transition from the non alternating to the dimerized 
regime can be followed by considering the second deriva- 
tive of the total energy (per carbon atom) with respect 
to S, E"(5) = E" a (S) + E'\(S). In the non alternating 
situation we have a minimum at <5 = 0, i.e. E"(0) > 0, 
while in the dimerized one we have two symmetrical min- 
ima at <5 ^ 0, i.e. <5 = is now a local maximum and 
E" (0) < 0. Therefore we can study the sign of the second 
derivative of the energy per carbon atom as a function 
of the number of atoms N: at the transition point E n 



5 



must be vanishing. In Table [ITj we report the values 
of EP\ per carbon atom for N = 10, 14, . . .30 together 
with E" „ per carbon atom (independent from TV); it is 
easily seen that the transition occurs between N = 14 
and N = 18. From these data it can be estimated that 
in order to have the transition at N = 30, we need to 
alter the ratio E" rT /E" 7V by a factor « 2, which cannot 
be considered a small adjustment of the model. 

In order to check the influence of geometry on these 
results, we choose the special case N = 18, and we com- 
pare energies obtained with the three different geometries 
of Fig. 0-|. The energies are shown in Table ttVl The 
geometry does not affect the dimerization in a sensible 
way, although the value of 5 is slightly reduced by going 
from (I) to (III). The values are the following: 



R0 = 1 .40 angstroms m = 51 2 



(I) ring: 

(II) cylinder: 

(III) annulene: 



Ro = 1.4067(6)1, 
R Q = 1.4064(7)1, 
R = 1.4063(1)1, 



6 = 0.0290(0)A 
S = 0.0276(1)1 
5 = 0.0271(6)1 



Table IV shows that the annulene geometry is variation- 
ally preferred. However, the a energy is completely unaf- 
fected by these geometry changes. The small differences 
in the values of Rq and S are only due to the 7r part of 
the Hamiltonian. 
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FIG. 7. Comparison of energies versus 8 at fixed 
Ro — 1AA. The symbols denote n energies per carbon atom 
for TV between 10 and 30, the solid line denotes the opposite 
of the a energy per carbon atom. Results are obtained with 
DMRG(2). 



TABLE III. Second derivatives of a and -k energies per 
carbon atom at zero dimerization, as obtained from DMRG(2) 
data. E'a does not depend on N. Derivatives are in atomic 
units. 



N 



K{8 = 0) 



K{8 = 0) 



10 
14 
IS 
22 
26 
30 



0.6082 



-0.44(0) 
-0.59(1) 
-0.74(9) 
-0.91(1) 
-1.06(7) 
-1.21(2) 



TABLE IV. Energy results for N = 18: comparison of 
7r-energies from DMRG(2) with different geometries, shown 
in Fig. 1, 2 and 3. Ro and 8 are in A, energies in eV. 



Ro 


8 


ring (I) 


cylinder (II) 


annulene (III) 


1.400 


0.020 


-44.826681 


-44.913191 


-44.951474 




0.025 


-44.952568 


-45.037015 


-45.074791 




0.030 


-45.098006 


-45.180696 


-45.217941 


1.405 


0.020 


-44.318921 


-44.404909 


-44.442656 




0.025 


-44.444324 


-44.528098 


-44.565347 




0.030 


-44.589043 


-44.670964 


-44.707697 


1.410 


0.020 


-43.815950 


-43.901321 


-43.938530 




0.025 


-43.940718 


-44.023869 


-44.060595 




0.030 


-44.084491 


-44.165918 


-44.202143 
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